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Wc simulate dopant profiles for phosphorus implantation into silicon using a new model 
for electronic stopping power. In this model, the electronic stopping power is factorized 
into a globally averaged effective charge Z^, and a local charge density dependent elec- 
tronic stopping power for a proton. There is only a single adjustable parameter in the 
model, namely the one electron radius r® which controls Z^. By fine tuning this pa- 
rameter, we obtain excellent agreement between simulated dopant profiles and the SIMS 
data over a wide range of energies for the channeling case. Our work provides a further 
example of implant species, in addition to boron and arsenic, to verify the validity of the 
electronic stopping power model and to illustrate its generality for studies of physical 
processes involving electronic stopping. 

1. Introduction 

Monte Carlo and molecular dynamics simulations of ion trajectories in a target ma- 
terial require a good description of the physics of electronic stopping in the high 
energy regime. The issue of the electronic stopping power is especially important 
when the target is crystalline and ions can propagate along preferred channel direc- 
tions in the lattice. For this case, electronic stopping becomes a dominant factor in 
determining final stopping ranges of the channeling ions. As is well known, the clas- 
sic Lindhard-Scharff-Schiott theoryS is applicable only to amorphous materials and 
underestimates electronic stopping along channeling directions. It thus tends to give 
an excessively high estimate of the energy threshold below which electronic stop- 
ping effects can safely be neglected when modeling ion implantation into crystalline 
materials. For the channeling case, a good understanding of electronic stopping is 
essential because the distribution of stopping ranges of ions can still be significantly 
affected by the electronic stopping power even at relatively low energies. 
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Recently we proposed a model for electronic stopping power for ion implanta- 
tion modcling.0 Based on the spirit of the Brandt-Kitagawa (BK) effective charge 
theory,0 it models the electronic stopping power for an ion in terms of two fac- 
tors: (i) a globally averaged effective charge taking into account effects of close and 
distant collisions by target electrons with the ion, and (ii) a local charge density 
dependent electronic stopping power for a proton. This model was implemented 
into both molecular dynamics and Monte Carlo simulations. The results of the 
dopant profile simulation for both boron and arsenic implants into crystalline and 
amorphous silicon have demonstrated that our model can successfully capture the 
physics of electronic stopping in ion implantation over a wide range of energies. 
The model is phenomenologically economical, i.e., it has only one tuning parame- 
ter, namely an averaged one electron radius r° which controls the effective charge 
of the ion. A single numerical value of this parameter was used for the simula- 
tion of both boron and arsenic implantation. Good agreement of dopant profiles 
with experimental profiles measured by secondary-ion mass spectroscopy (SIMS) 
was achieved for both species with this single r® numerical value. 

We note that the BK theory uses a statistical model for the partially ionized 
projectile and does not account for shell structure. Therefore, it can only provide 
an averaged description of electronic stopping power as a function of the projectile 
atomic number Z\ .□ The experimentally observed Z\ oscillations in electronic stop- 
ping are a complex phenomenon, attributable JjQ-thc electronic shell structures of 
both the incident ion and the target atom.B&LBlj On account of the dependence 
of the Z\ oscillation on both the ion and target material, we expect that the pa- 
rameter r® can be tuned to different numerical values for different combinations 
of implant species and substrate material. This fine tuning can be viewed as a 
phenomcnological procedure to incorporate the physics of Z\ oscillations. In the 
present work, we will verify this phenomenological approach for phosphorus im- 
plantation into silicon and show that our electronic stopping power can successfully 
model channeling of phosphorus implants into single-crystal silicon, thus extending 
our electronic stopping power model to the case of phosphorus-on-silicon implants. 
This further illustrates the potential wide applicability of the model in studies of 
physical processes that involve electronic stopping. 

The paper is organized as follows. In Sec. 2 we summarize the main features of 
our electronic stopping power model and briefly compare it with other models that 
have been used in various Monte Carlo simulations. Atomic units e = h = m e = 1 
are used in the description of our model. In_Sec. 3 we describe the implementation 
of our model into a MARLOWE platform,t2l and into a molecular dynamics (MD) 
based implant simulator.til We then present a comparison between our simulation 
results and SIMS data for phosphorus implantation into silicon. In Sec. 4 we make 
concluding remarks. 
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2. The Model 



In our model Ja based on an effective charge scaling argument, the electronic stopping 
power of an ion can be factorized into two components. One is the effective charge 
Z* of the ion of velocity v\ , which can be expressed as 



Z{ = Zi7(vi,r°), 



(1) 



where Z\ is the atomic number of the ion and 7(1*1, r°) is the fractional effec- 
tive charge of the ion. The second is the charge density dependent electronic 
stopping power for a proton S p (vi,r s ). Here r s is the one electron radius, r s = 
[3/ (47T/9(x))] 1 ^ 3 , where p(x) is the charge density of the target. In our treatment 
7(1*1,7*;;) does not depend on the local charge density, instead, it is controlled by 
the parameter r®, which is the only adjustable parameter in the model. 

After taking account of the energy loss of the ion in soft, distant collisions with 
target electrons and the energy loss in hard close collisions, the BK analysis gives 
the fractional effective charge 



1 (v 1 ,r s ) = q(v 1 ,r s ) + C[l-q(v u r° s )}ln 



(2) 



Here C is weakly dependent on the target and has a numerical value of about 1/2. 
Below, it is set to be one-half. The ion size parameter A(i*i,rg) is 



2ao[l-g(t>i,r°)] 2 /3 



zni-(i-<?M))/7] 



(3) 



which is used in the statistical model to describe the partially ionized projectile. 
Here ao — 0.24 005 and the ionization fraction q(vi,r®) obeys the following scaling 

<?(«!, r°) = l-exp[-0.95(y r (vi,r°)-0.07)]. (4) 

This scaling was condensed from extensive experimental data for ions 3 < Z 1 < 920 
The reduced relative velocity y r (vi,r®) is defined as 



Vr(vi,r° s ) 



v B Z 1 



2/3 



(5) 



where vg is the Bohr velocity. Underlying Eq. (||) is the stripping criterion that the 
electrons of the ion which have an orbital velocity lower than the relative velocity 
between the ion and the electrons in the medium_are stripped off. Averaging relative 
velocity over the conduction electrons leads to:ll3 



v r (vi,r° s ) 



1*1 I 1 — I l,,r : 



I On 3l*F ( 2vf 



154 



for vi < vf 



(6) 
(7) 
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In our modclja the Fermi velocity is related to r®: vf — l/far®), a = [4/(9 7r )] 1// ' 3 - 

The electronic stopping power for a proton in our model uses results derived 
from a nonlinear density-functional formalismiia 



.S„(r l . ; .. = -('gj a G(,.!. (M 



where 

1 



dx J R 37T 



In 1 



(9) 



1 + ar s /ir 

is the Ritchie formula for the energy loss per unit path length of a proton moving 
at velocity V\ in the electron gas derived from a linear response theory. The correc- 
tion factor G{r s ) is a computationally convenient way to incorporate the density- 
functional results and it has the following form 

G{r s ) = 1.00 + 0.717r s - 0.125^ - 0.0124r 3 + 0.00212^ (10) 

for r s < 6. It should be emphasized that S p (vi,r s ) in our model depends on the 
local charge density, p(x) = 3/(47rr s (x) 3 ). We note in passing that the density- 
functional result gives a better estimate than the linear response (dielectric) result 



for t 
data. 



le stopping powers as demonstrated in the comparisons with experimental 
^0 The charge density, p(x), for atoms in crystal silicon in our model uses 



the solid-state Hartree-Fock atomic charge distribution.0 In this approximation, 
there is about a one electron charge not accounted for inside the spherical muffin- 
tin. This small amount of charge is distributed between the maximal collision 
distance used in our Monte Carlo simulations and the muffin-tin radius J§ within the 
MD scheme it provides a background charge experienced by the ion when further 
than the muffin-tin radius from any silicon atomJl^l 

In our simulation, the electronic stopping power is evaluated continuously along 
the path the ion traverses through regions of varying charge density. The energy 
loss is computed by 

AE e = [ [Z l7 (v lir ° s )] 2 S p {v u r s {x))dx. (11) 

J ion path 

Finally, we discuss the difference between our approach and other electronic 
stopping power models used in Monte Carlo simulations based on the MARLOWE 
platform. There is a purely nonlocal version of the BK theory implemented into 
MARLOWE J13 where the effective charge and the stopping power for a proton both 
depend on a single nonlocal parameter, i.e., the averaged one electron radius. The 
energy loss for well-channeled ions in the keV region in this approach indicated that 
a correc t d e n sity distribution is necessary to model the electronic stopping in the 
channel. tllta Later, an implementation of a purely local version of the BK theory to 
model boron implants into (100) single-crystal silicon produced dopant profiles in 
good agreement with the SIMS data.E3 It showed a marked improvement in modeling 
dopant profiles in the channel over other electronic stopping power models, such as 



Simulation of Phosphorus Implantation Into Silicon 5 



Lindhard and Scharff J23 FirsovJ^ and the above nonlocal implementation. However, 
this purely local version was unable to model the electronic storming either for boron 
implants into the (110) axial channel, or for arsenic implants.E3 

In our model, the effective charge is a nonlocal quantity, neither explicitly depen- 
dent on the impact parameter nor on the charge distribution, while the stopping 
power for a proton depends on the local charge density of silicon. It has been 
shown that it can successfully model both boron and arsenic implants into silicon 
over a wide range of energies and in different channeling and off-axis directions of 
incidencej§ : In the following section, we will demonstrate that this model can be 
further extended to phosphorus implantation into silicon. In light of Z\ oscillations, 
r° will take a slightly different numerical value from the one used in the simulations 
for boron and arsenic implantation.B 

3. Monte Carlo and Molecular Dynamics Simulation Results 

We have implemented the electronic stopping power model into the Monte Carlo 
simulation platform MARLOWE, which utilizes the binary collision approximation 
(BCA) to simulate the trajectories of energetic ions in crystalline or amorphous 
materials. EH For the present work, we used an extended version of MARLOWE 
with enhanced, capabilities for modeling ion implant into silicon. This code, UT- 
MARLOWEjlJ was selected for its more versatile features, including variance re- 
duction for more efficient calculation, and the incorporation of important implant 
parameters, e.g., tilt and rotation angles, the thickness of the native oxide layers, 
beam divergence, and wafer temperature. The electronic stopping power model is 
also incorporated into an MD based implant simulator, REEDO This allows us 
to calculate profiles using a more realistic description of atomic collisions during 
implants than that provided by the BCA. It also gives us an additional verification 
of the model, and shows that it is independent of the simulation platform. 

For the purpose of verifying the electronic stopping model, the option of the 
cumulative damage model in the UT-MARLOWE code was turned off in our sim- 
ulations. This is a phenomenological model for estimating defect production and 
recombination rates. The REED cumulative damage model was also disabled. In 
order to minimize cumulative damage effects in the dopant profiles, simulations 
were performed for low to moderate dose (10 13 cm -2 to 3 x 10 13 cm -2 ) phospho- 
rus implants. Individual ion trajectories were simulated and the overlapping of the 
damage caused by different individual ion cascades was neglected. Also, use was 
made of 1° initial beam divergence, a 16 A native oxide surface layer, and 300 K 
wafer temperature. These parameters were not cited in the experimental references 
but are believed to be typical for silicon doping implants. 

The best suited data for verifying the electronic stopping model are the on-axis 
or (100) channeling implants. These data are highly sensitive to electronic stopping 
and are relatively insensitive to other effects, especially at energies of 100 keV and 
above where electronic stopping becomes the dominant energy loss mechanism. The 
free parameter in the model r® was adjusted to 1.148 A to yield the best results in 



6 D. Cat, C. M. Snell, K. M. Beardmore & N. Gr0nbech- Jensen 



overall comparison between the BCA and SIMS profiles. A fixed numerical value 
of r® was employed for all energies and directions of incidence. As commented 
beforep this value is physically reasonable for silicon. We note that there is only 
a 3.5% difference between this value of r|? and the value of 1.109 A previously 
used to obtain the best overall agreement with experimental data for boron and 
arsenic.cl This difference is quite small, and our result indicates that, in general, 
Z\ oscillations can be accounted for by fine-tuning of r° for a given combination 
of implant species and material. A slightly larger value of 1.217 A was used for r° 
within the MD simulations. This difference reflects the fact that the description of 
ion channelling and energy loss differs between the BCA and MD schemes, and that 
while fitting the electronic stopping model we are also compensating for deficiencies 
in the implant simulation models. 



i r 




1000 1500 
Depth (A) 

Fig. 1. Calculated and experimental dopant profiles due to 15 keV phosphorus implant for the 
(100) direction with zero tilt and rotation angles; the dose is 10 13 cm -2 . 




Fig. 2. 
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Calculated and experimental dopant profiles due to 50 keV phosphorus implant for the 



(100) direction with zero tilt and rotation angles; the dose is 10 cm - the SIMS profile for a 
dose of 3 X 10 13 cm -2 is also shown to illustrate damage effects. 
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We now turn to our Monte Carlo and MD dopant profile simulation results. For 
comparison we have digitized the SIMS data for phosphorus reported in Refs. 9 
and 25. The simulated phosphorus dopant profiles are shown in comparison to the 
SIMS data for the channeling case in Figs. ^ |[ and |]. Here, the implant energy 
ranges from 15 keV to 200 keV. The incidence is along the (100) direction with 
zero tilt and zero rotation angles. The dose is 10 13 cm" 2 or 3 x 10 13 cm" 2 . The 
simulations show a close fit to the dopant distributions with depth and describe 
especially well the slope of the channeling tails of the dopant profiles. Some un- 
certainty is introduced into these comparisons by the fact that the exact implant 
parameters employed in the experiments are not known. Ion channeling ranges for 
on-axis implants may be sensitive to variations in these parameters, especially to 
the beam divergence. Fig. || illustrates the effect of changing the assumed beam 
divergence from 1° to 0°. As expected, the lower divergence gives a shallower slope 
and deeper penetration in the channeling tail. Best agreement with the SIMS profile 
is achieved for a divergence angle of 1°, which is quite reasonable for commercial 
implant machines. This value was employed in all other calculations presented here. 

Another source of uncertainty is the possible influence of disorder introduced by 
damage to the silicon, which can reduce the average range of channeled ions. This 
effect was ignored in our calculations. The importance of damage can be roughly 
evaluated by comparing the measured SIMS profiles at different doses, as in Figs. || 
and |. It is seen that the profiles at the higher dose are similar the lower dose 
profiles, but that a three-fold increase in dose results in an almost five-fold increase 
in dopant concentration at the peak while only doubling the concentration in the 
tail. However the depth of the peak and end-of-range of the profile are unaltered 
by the addition of damage. Hence, we can compare the 'zero damage' simulations 
to the higher dose (3 x 10 13 cm" 2 ) SIMS data, but must be aware that the slope 
of the channelling tail may not match - the MD profile in Fig. [| shows this effect. 
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Fig. 4. Calculated and experimental dopant profiles due to 200 keV phosphorus implant for the 
(100) direction with zero tilt and rotation angles; the dose is 3 X 10 13 cm -2 . 




2000 4000 6000 8000 10000 
Depth (A) 

Fig. 5. Comparison of profiles calculated with 0° and 1° beam divergence; other parameters are 
as Fig. | 

We conclude that the overall agreement between the simulations and experimental 
data for the channeling cases is excellent over a wide range of energies - especially 
considering that the electronic stopping model implemented in the BCA platform 
has been tuned by only 3.5% in the parameter, ri°K 

Figs. H and |t] show the simulated dopant profiles in the off-axis directions, 
10° tilt and 15° rotation for 100 keV, and 8° tilt and 18° rotation for 200 keV, 
respectively. The BCA calculated profiles show less good agreement with SIMS 
data than for channeling cases, with systematically deeper penetration and the 
concentration peak shifted by about 25% - a similar, but smaller shift is also seen 
for the 200 keV on-axis implant (Fig. [|). 

The explanation for this discrepancy lies in the different roles of energy loss and 
scattering mechanisms for the off-axis implants. The average final stopping range 
for ions implanted in off-axis directions is controlled primarily by nuclear scattering, 
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Fig. 6. Calculated and experimental dopant profiles due to 100 keV phosphorus implant for the 
(100) direction with 10° tilt and 15° rotation; the dose is 10 13 cm~ 2 . 




4000 6000 8000 10000 
Depth (A) 

Fig. 7. Calculated and experimental dopant profiles due to 200 keV phosphorus implant for the 
(100) direction with 8° tilt and 18° rotation; the dose is 3 X 10 13 cm -2 . 

and also by inelastic energy loss0 during ion-atom interactions; it is less sensitive to 
the electronic stopping model (except at very high energies above those investigated 
in this study). Accurate prediction of the off-axis profiles thus requires an accurate 
model for the atom-specific interatomic potential of the two species involved in 
the implant, plus a description of inelastic energy loss. We do not currently have 
an atom-specific potential model for phosphorus and silicon in the BCA platform, 
and there is no attempt to account for inelastic energy loss due to momentum 
transfer between electrons during collisions. The significance of these omissions 
in the BCA model was examined by performing separate calculations for random 
implant directions and for an implant into amorphous silicon (not shown here). 
These cases are completely dominated by nuclear scattering and are insensitive to 
electronic stopping in this energy range. The penetration ranges were too large by 
about 20%, verifying that nuclear scattering and inelastic energy loss are largely 
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responsible for the over prediction of ranges seen in the off- axis calculations. 

The MD platform (REED) uses pair-specific potentialsEl for all ion-silicon in- 
teractions, and includes an inelastic energy loss model.E] The off-axis profiles cal- 
culated by MD have the correct peak position and shallower penetration than the 
BCA profiles, and provide a better match to the SIMS data. The BCA predictions 
could probably be improved by the relatively simple addition of the two models 
mentioned above to the BCA platform, although that is outside the scope of this 
paper. Profiles calculated by both simulation models show an increased amount of 
channelling relative to the SIMS profiles. The channelling tail is due to incident ions 
being scattered into channelling directions in near-surface collisions. The number 
hard collisions near the surface due to off-axis implants will produce a significant 
amount of damage, even at a dose as low as I0 13 cm -2 . Hence we should expect 
the 'zero damage' calculated profiles to show exaggerated channelling relative to 
the experimental data. 

4. Conclusion 

We have used our electronic stopping power model to simulate dopant profiles for 
phosphorus implantation into silicon. To account for Z\ oscillations, we have ap- 
propriately fine-tuned the single adjustable parameter r® in our model to match the 
phosphorus-silicon case. The numerical value of r® is slightly greater than the one 
used for boron and arsenic simulations. Using this r° = 1.148 A, our BCA results 
show excellent agreement between the simulated dopant profiles and the SIMS data 
over a wide range of energies for the channeling case. Less accurate but satisfac- 
tory results are obtained for off-axis implants. Detailed agreement in the off-axis 
direction would require additional model development for ion-silicon interactions. 
We have also implemented the stopping model in an MD based simulator, using 
a slightly larger value for r®. Using the MD implant model we achieve excellent 
agreement with SIMS data for both on-axis, and off-axis implants. In summary, 
we have successfully extended our electronic stopping power model to encompass 
phosphorus implantation into crystalline silicon. We have also indicated how to 
incorporate Z\ oscillations with a simple phenomenological approach. We have 
provided a further example of implant species to verify validity of the model and 
to demonstrate its generality for studies of physical processes involving electronic 
stopping. 
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